program PS

        ! "use" is the inheritance mechanism for F90:
        use mConstants
        use mObject
        use mCosmology
        use mLPS
        use mPSFittingFunctions
        use mMisc

        implicit none

        real(D) :: smithp3d, linearp3d
        real(D) :: lnk, dlnk, k, x, w, w2, p, dsig2, sig2, twopisq
        integer :: i

        ! om_m, om_lam, om_b, h0, As, n, w, wa
        call def_universe(0.27_D,0.73_D,0.046_D,0.72_D,0.81_D,1.0_D,-1._D,0._D)

        lnk = -7.
        dlnk = 0.01
        sig2 = 0.
        twopisq = 19.7392088022

        write(6,*) 'k(h/Mpc)  W_sigma  Delta3D  d(sigma_8^2)'
        do i = 1, 1000

          k = exp(lnk)
          x = k*8./0.72_D  ! k in h/Mpc, 8 in h^-1 Mpc, product dimensionless
          w = (sin(x)-x*(cos(x))) * 3 / (x**3)
          w2 = w*w;
          !p = smithp3d(k,0.) ! = P(k)*k^3/twopisq
          p = linearp3d(k,0.) * k**3 / twopisq

          dsig2 = w2*p*dlnk

          write(6,'(4e12.3)') k,w,p,dsig2

          sig2 = sig2 + dsig2
          lnk = lnk + dlnk
        end do

        write(*,*) 'sigma_8^2 = ',sig2
        write(*,*) 'sigma_8 = ',sqrt(sig2)

end program PS
